
##########################################################
## SETUP (Load required packages, set working directory)
##########################################################

library(survival)
library(coxed)
library(dplyr)
library(tidyr)
library(ggplot2)
library(foreign)

setwd("~/Dropbox/Projekt/Corona/JPP/JPP Replication/Dataverse")

##########################################################
## MODELS A & A+ (Baseline variables)
##########################################################

## A

JPP_covid19_2022_CoxED_A <- read.csv("JPP_covid19_2022_CoxED_A.csv")

covid.surv_A <- Surv(time = JPP_covid19_2022_CoxED_A$date1, time2 = JPP_covid19_2022_CoxED_A$date2, event = JPP_covid19_2022_CoxED_A$schoolsclosed)
covid.cox_A <- coxph(covid.surv_A ~ v2x_polyarchy2019 + wbgi_gee2019 + ln_cases1 + reg_any_case + weekend + strata(date_firstcase) + as.factor(region) + cluster(ccode) , method = "breslow", data = JPP_covid19_2022_CoxED_A)
summary(covid.cox_A)

edA <- coxed(covid.cox_A, method = "npsf",  id = JPP_covid19_2022_CoxED_A$ccode, data=JPP_covid19_2022_CoxED_A)
summary(edA, stat="mean")

# A-D:
me_A_D <- coxed(covid.cox_A, method = "npsf",  id = JPP_covid19_2022_CoxED_A$ccode, bootstrap = TRUE, B = 200, maxit = 30,
                newdata = dplyr::mutate(JPP_covid19_2022_CoxED_A, v2x_polyarchy2019 = JPP_covid19_2022_CoxED_A$lowD),
                newdata2 = dplyr::mutate(JPP_covid19_2022_CoxED_A, v2x_polyarchy2019 = JPP_covid19_2022_CoxED_A$highD))
summary(me_A_D, stat="mean")

# A-SC: 
me_A_SC <- coxed(covid.cox_A, method = "npsf",  id = JPP_covid19_2022_CoxED_A$ccode, bootstrap = TRUE, B = 200, maxit = 30,
                 newdata = dplyr::mutate(JPP_covid19_2022_CoxED_A, wbgi_gee2019   = lowSC),
                 newdata2 = dplyr::mutate(JPP_covid19_2022_CoxED_A, wbgi_gee2019   = highSC))
summary(me_A_SC, stat="mean")

## A + Interaction

covid.cox_A.int <- coxph(covid.surv_A ~ v2x_polyarchy2019 * wbgi_gee2019 + ln_cases1 + reg_any_case + weekend + strata(date_firstcase) + as.factor(region) + cluster(ccode) , method = "breslow", data = JPP_covid19_2022_CoxED_A)
summary(covid.cox_A.int)

# A-D-low-SC:
me_A_int_DlowSC <- coxed(covid.cox_A.int, method = "npsf",  id = JPP_covid19_2022_CoxED_A$ccode, bootstrap = TRUE, B = 200, maxit = 30,
                         newdata = dplyr::mutate(JPP_covid19_2022_CoxED_A, v2x_polyarchy2019 = JPP_covid19_2022_CoxED_A$lowD, wbgi_gee2019 = JPP_covid19_2022_CoxED_A$lowSC),
                         newdata2 = dplyr::mutate(JPP_covid19_2022_CoxED_A, v2x_polyarchy2019 = JPP_covid19_2022_CoxED_A$highD, wbgi_gee2019 = JPP_covid19_2022_CoxED_A$lowSC))
summary(me_A_int_DlowSC, stat="mean")

#A-D-high-SC:
me_A_int_DhighSC <- coxed(covid.cox_A.int, method = "npsf",  id = JPP_covid19_2022_CoxED_A$ccode, bootstrap = TRUE, B = 200, maxit = 30,
                          newdata = dplyr::mutate(JPP_covid19_2022_CoxED_A, v2x_polyarchy2019 = JPP_covid19_2022_CoxED_A$lowD, wbgi_gee2019 = JPP_covid19_2022_CoxED_A$highSC),
                          newdata2 = dplyr::mutate(JPP_covid19_2022_CoxED_A, v2x_polyarchy2019 = JPP_covid19_2022_CoxED_A$highD, wbgi_gee2019 = JPP_covid19_2022_CoxED_A$highSC))
summary(me_A_int_DhighSC, stat="mean")

# A-SC-low-D:
me_A_int_DlowSC <- coxed(covid.cox_A.int, method = "npsf",  id = JPP_covid19_2022_CoxED_A$ccode, bootstrap = TRUE, B = 200, maxit = 30,
                         newdata = dplyr::mutate(JPP_covid19_2022_CoxED_A, v2x_polyarchy2019 = JPP_covid19_2022_CoxED_A$lowD, wbgi_gee2019 = JPP_covid19_2022_CoxED_A$lowSC),
                         newdata2 = dplyr::mutate(JPP_covid19_2022_CoxED_A, v2x_polyarchy2019 = JPP_covid19_2022_CoxED_A$lowD, wbgi_gee2019 = JPP_covid19_2022_CoxED_A$highSC))
summary(me_A_int_DlowSC, stat="mean")

#A-SC-high-D:
me_A_int_DhighSC <- coxed(covid.cox_A.int, method = "npsf",  id = JPP_covid19_2022_CoxED_A$ccode, bootstrap = TRUE, B = 200, maxit = 30,
                          newdata = dplyr::mutate(JPP_covid19_2022_CoxED_A, v2x_polyarchy2019 = JPP_covid19_2022_CoxED_A$highD, wbgi_gee2019 = JPP_covid19_2022_CoxED_A$lowSC),
                          newdata2 = dplyr::mutate(JPP_covid19_2022_CoxED_A, v2x_polyarchy2019 = JPP_covid19_2022_CoxED_A$highD, wbgi_gee2019 = JPP_covid19_2022_CoxED_A$highSC))
summary(me_A_int_DhighSC, stat="mean")


##########################################################
## MODELS B & B+ (Baseline variables)
##########################################################

## B (A + Additional controls)

JPP_covid19_2022_CoxED_B <- read.csv("JPP_covid19_2022_CoxED_B.csv")

covid.surv_B <- Surv(time = JPP_covid19_2022_CoxED_B$date1, time2 = JPP_covid19_2022_CoxED_B$date2, event = JPP_covid19_2022_CoxED_B$schoolsclosed)
covid.cox_B <- coxph(covid.surv_B ~ v2x_polyarchy2019 + wbgi_gee2019 + ln_cases1 + reg_any_case + weekend + ln_gdpcap + wdi_beds_per1000c + wdi_pop14 + wdi_popurb + acd_res_nuclstem + strata(date_firstcase) + as.factor(region) + cluster(ccode) , method = "breslow", data = JPP_covid19_2022_CoxED_B)
summary(covid.cox_B)

ed_B <- coxed(covid.cox_B, method = "npsf",  id = JPP_covid19_2022_CoxED_B$ccode, data=JPP_covid19_2022_CoxED_B)
summary(ed_B, stat="mean")

# B-D:
me_B_D <- coxed(covid.cox_B, method = "npsf",  id = JPP_covid19_2022_CoxED_B$ccode, bootstrap = TRUE, B = 200, maxit = 30,
                newdata = dplyr::mutate(JPP_covid19_2022_CoxED_B, v2x_polyarchy2019 = JPP_covid19_2022_CoxED_B$lowD),
                newdata2 = dplyr::mutate(JPP_covid19_2022_CoxED_B, v2x_polyarchy2019 = JPP_covid19_2022_CoxED_B$highD))
summary(me_B_D, stat="mean")

# B-SC: 
me_B_SC <- coxed(covid.cox_B, method = "npsf",  id = JPP_covid19_2022_CoxED_B$ccode, bootstrap = TRUE, B = 200, maxit = 30,
                 newdata = dplyr::mutate(JPP_covid19_2022_CoxED_B, wbgi_gee2019 = lowSC),
                 newdata2 = dplyr::mutate(JPP_covid19_2022_CoxED_B, wbgi_gee2019 = highSC))
summary(me_B_SC, stat="mean")


## B + Interaction

covid.cox_B.int <- coxph(covid.surv_B ~ v2x_polyarchy2019 * wbgi_gee2019 + ln_cases1 + reg_any_case + weekend + ln_gdpcap + wdi_beds_per1000c + wdi_pop14 + wdi_popurb + acd_res_nuclstem + strata(date_firstcase) + as.factor(region) + cluster(ccode) , method = "breslow", data = JPP_covid19_2022_CoxED_B)
summary(covid.cox_B.int)

# B-D-low-SC:
me_B_int_DlowSC <- coxed(covid.cox_B.int, method = "npsf",  id = JPP_covid19_2022_CoxED_B$ccode, bootstrap = TRUE, B = 200, maxit = 30,
                         newdata = dplyr::mutate(JPP_covid19_2022_CoxED_B, v2x_polyarchy2019 = JPP_covid19_2022_CoxED_B$lowD, wbgi_gee2019 = JPP_covid19_2022_CoxED_B$lowSC),
                         newdata2 = dplyr::mutate(JPP_covid19_2022_CoxED_B, v2x_polyarchy2019 = JPP_covid19_2022_CoxED_B$highD, wbgi_gee2019 = JPP_covid19_2022_CoxED_B$lowSC))
summary(me_B_int_DlowSC, stat="mean")

#B-D-high-SC:
me_B_int_DhighSC <- coxed(covid.cox_B.int, method = "npsf",  id = JPP_covid19_2022_CoxED_B$ccode, bootstrap = TRUE, B = 200, maxit = 30,
                          newdata = dplyr::mutate(JPP_covid19_2022_CoxED_B, v2x_polyarchy2019 = JPP_covid19_2022_CoxED_B$lowD, wbgi_gee2019 = JPP_covid19_2022_CoxED_B$highSC),
                          newdata2 = dplyr::mutate(JPP_covid19_2022_CoxED_B, v2x_polyarchy2019 = JPP_covid19_2022_CoxED_B$highD, wbgi_gee2019 = JPP_covid19_2022_CoxED_B$highSC))
summary(me_B_int_DhighSC, stat="mean")

# B-SC-low-D:
me_B_int_DlowSC <- coxed(covid.cox_B.int, method = "npsf",  id = JPP_covid19_2022_CoxED_B$ccode, bootstrap = TRUE, B = 200, maxit = 30,
                         newdata = dplyr::mutate(JPP_covid19_2022_CoxED_B, v2x_polyarchy2019 = JPP_covid19_2022_CoxED_B$lowD, wbgi_gee2019   = JPP_covid19_2022_CoxED_B$lowSC),
                         newdata2 = dplyr::mutate(JPP_covid19_2022_CoxED_B, v2x_polyarchy2019 = JPP_covid19_2022_CoxED_B$lowD, wbgi_gee2019   = JPP_covid19_2022_CoxED_B$highSC))
summary(me_B_int_DlowSC, stat="mean")

#B-SC-high-D:
me_B_int_DhighSC <- coxed(covid.cox_B.int, method = "npsf",  id = JPP_covid19_2022_CoxED_B$ccode, bootstrap = TRUE, B = 200, maxit = 30,
                          newdata = dplyr::mutate(JPP_covid19_2022_CoxED_B, v2x_polyarchy2019 = JPP_covid19_2022_CoxED_B$highD, wbgi_gee2019   = JPP_covid19_2022_CoxED_B$lowSC),
                          newdata2 = dplyr::mutate(JPP_covid19_2022_CoxED_B, v2x_polyarchy2019 = JPP_covid19_2022_CoxED_B$highD, wbgi_gee2019   = JPP_covid19_2022_CoxED_B$highSC))
summary(me_B_int_DhighSC, stat="mean")


##########################################################
## MODELS C & D (Vanhanen: Competition and Participation)
##########################################################

## C

JPP_covid19_2022_CoxED_C <- read.csv("JPP_covid19_2022_CoxED_C.csv")

covid.surv_C <- Surv(time = JPP_covid19_2022_CoxED_C$date1, time2 = JPP_covid19_2022_CoxED_C$date2, event = JPP_covid19_2022_CoxED_C$schoolsclosed)
covid.cox_C <- coxph(covid.surv_C ~ van_comp2018 + van_part2018 + wbgi_gee2019 + ln_cases1 + reg_any_case + weekend + strata(date_firstcase) + as.factor(region) + cluster(ccode) , method = "breslow", data = JPP_covid19_2022_CoxED_C)
summary(covid.cox_C)

ed_C <- coxed(covid.cox_C, method = "npsf",  id = JPP_covid19_2022_CoxED_C$ccode, data=JPP_covid19_2022_CoxED_C)
summary(ed_C, stat="mean")

# C: Compare Low P and Low C to Low C and High P
me_C_auth_Cl <- coxed(covid.cox_C, method = "npsf",  id = JPP_covid19_2022_CoxED_C$ccode, bootstrap = TRUE, B = 200, maxit = 30,
                      newdata = dplyr::mutate(JPP_covid19_2022_CoxED_C, van_comp2018 = JPP_covid19_2022_CoxED_C$lowC, van_part2018 = JPP_covid19_2022_CoxED_C$lowP),
                      newdata2 = dplyr::mutate(JPP_covid19_2022_CoxED_C, van_comp2018 = JPP_covid19_2022_CoxED_C$lowC, van_part2018 = JPP_covid19_2022_CoxED_C$highP))
summary(me_C_auth_Cl, stat="mean")

# C: Compare Low P and Low C to High C and High P
me_C_auth_dem <- coxed(covid.cox_C, method = "npsf",  id = JPP_covid19_2022_CoxED_C$ccode, bootstrap = TRUE, B = 200, maxit = 30,
                       newdata = dplyr::mutate(JPP_covid19_2022_CoxED_C, van_comp2018 = JPP_covid19_2022_CoxED_C$lowC, van_part2018 = JPP_covid19_2022_CoxED_C$lowP),
                       newdata2 = dplyr::mutate(JPP_covid19_2022_CoxED_C, van_comp2018 = JPP_covid19_2022_CoxED_C$highC, van_part2018 = JPP_covid19_2022_CoxED_C$highP))
summary(me_C_auth_dem, stat="mean")

# C: Compare High P and Low C to High C and High P (not reported)
me_C_Cl_dem <- coxed(covid.cox_C, method = "npsf",  id = JPP_covid19_2022_CoxED_C$ccode, bootstrap = TRUE, B = 200, maxit = 30,
                     newdata = dplyr::mutate(JPP_covid19_2022_CoxED_C, van_comp2018 = JPP_covid19_2022_CoxED_C$lowC, van_part2018 = JPP_covid19_2022_CoxED_C$highP),
                     newdata2 = dplyr::mutate(JPP_covid19_2022_CoxED_C, van_comp2018 = JPP_covid19_2022_CoxED_C$highC, van_part2018 = JPP_covid19_2022_CoxED_C$highP))
summary(me_C_Cl_dem, stat="mean")


## D (C + Additional controls)

JPP_covid19_2022_CoxED_D <- read.csv("JPP_covid19_2022_CoxED_D.csv")

covid.surv_D <- Surv(time = JPP_covid19_2022_CoxED_D$date1, time2 = JPP_covid19_2022_CoxED_D$date2, event = JPP_covid19_2022_CoxED_D$schoolsclosed)
covid.cox_D <- coxph(covid.surv_D ~ van_comp2018 + van_part2018 + wbgi_gee2019 + ln_cases1 + reg_any_case + weekend + ln_gdpcap + wdi_beds_per1000c + wdi_pop14 + wdi_popurb + acd_res_nuclstem + strata(date_firstcase) + as.factor(region) + cluster(ccode) , method = "breslow", data = JPP_covid19_2022_CoxED_D)
summary(covid.cox_D)

ed_D <- coxed(covid.cox_D, method = "npsf",  id = JPP_covid19_2022_CoxED_D$ccode, data=JPP_covid19_2022_CoxED_D)
summary(ed_D, stat="mean")

# D: Compare Low P and Low C to Low C and High P
me_D_auth_el <- coxed(covid.cox_D, method = "npsf",  id = JPP_covid19_2022_CoxED_D$ccode, bootstrap = TRUE, B = 200, maxit = 30,
                      newdata = dplyr::mutate(JPP_covid19_2022_CoxED_D, van_comp2018 = JPP_covid19_2022_CoxED_D$lowC, van_part2018 = JPP_covid19_2022_CoxED_D$lowP),
                      newdata2 = dplyr::mutate(JPP_covid19_2022_CoxED_D, van_comp2018 = JPP_covid19_2022_CoxED_D$lowC, van_part2018 = JPP_covid19_2022_CoxED_D$highP))
summary(me_D_auth_el, stat="mean")

# D: Compare Low P and Low C to High C and High P
me_D_auth_dem <- coxed(covid.cox_D, method = "npsf",  id = JPP_covid19_2022_CoxED_D$ccode, bootstrap = TRUE, B = 200, maxit = 30,
                       newdata = dplyr::mutate(JPP_covid19_2022_CoxED_D, van_comp2018 = JPP_covid19_2022_CoxED_D$lowC, van_part2018 = JPP_covid19_2022_CoxED_D$lowP),
                       newdata2 = dplyr::mutate(JPP_covid19_2022_CoxED_D, van_comp2018 = JPP_covid19_2022_CoxED_D$highC, van_part2018 = JPP_covid19_2022_CoxED_D$highP))
summary(me_D_auth_dem, stat="mean")

# D: Compare High P and Low C to High C and High P (not reported)
me_D_el_dem <- coxed(covid.cox_D, method = "npsf",  id = JPP_covid19_2022_CoxED_D$ccode, bootstrap = TRUE, B = 200, maxit = 30,
                     newdata = dplyr::mutate(JPP_covid19_2022_CoxED_D, van_comp2018 = JPP_covid19_2022_CoxED_D$lowC, van_part2018 = JPP_covid19_2022_CoxED_D$highP),
                     newdata2 = dplyr::mutate(JPP_covid19_2022_CoxED_D, van_comp2018 = JPP_covid19_2022_CoxED_D$highC, van_part2018 = JPP_covid19_2022_CoxED_D$highP))
summary(me_D_el_dem, stat="mean")


##########################################################
## MODELS E & F (fh_ipolity2 and icrg_qog)
##########################################################

## E

JPP_covid19_2022_CoxED_E <- read.csv("JPP_covid19_2022_CoxED_E.csv")

covid.surv_E <- Surv(time = JPP_covid19_2022_CoxED_E$date1, time2 = JPP_covid19_2022_CoxED_E$date2, event = JPP_covid19_2022_CoxED_E$schoolsclosed)
covid.cox_E <- coxph(covid.surv_E ~ fh_ipolity2 + icrg_qog +ln_cases1 + reg_any_case + weekend + strata(date_firstcase) + as.factor(region) + cluster(ccode) , method = "breslow", data = JPP_covid19_2022_CoxED_E)
summary(covid.cox_E)

edG <- coxed(covid.cox_E, method = "npsf",  id = JPP_covid19_2022_CoxED_E$ccode, data=JPP_covid19_2022_CoxED_E)
summary(edG, stat="mean")

# E-D:
me_E_D <- coxed(covid.cox_E, method = "npsf",  id = JPP_covid19_2022_CoxED_E$ccode, bootstrap = TRUE, B = 200, maxit = 30,
                newdata = dplyr::mutate(JPP_covid19_2022_CoxED_E, fh_ipolity2 = JPP_covid19_2022_CoxED_E$lowFPDem),
                newdata2 = dplyr::mutate(JPP_covid19_2022_CoxED_E, fh_ipolity2 = JPP_covid19_2022_CoxED_E$highFPDem))
summary(me_E_D, stat="mean")

# E-QOG: 
me_E_QOG <- coxed(covid.cox_E, method = "npsf",  id = JPP_covid19_2022_CoxED_E$ccode, bootstrap = TRUE, B = 200, maxit = 30,
                  newdata = dplyr::mutate(JPP_covid19_2022_CoxED_E, icrg_qog = lowQOG),
                  newdata2 = dplyr::mutate(JPP_covid19_2022_CoxED_E, icrg_qog = highQOG))
summary(me_E_QOG, stat="mean")


## E + Interaction

covid.cox_E.int <- coxph(covid.surv_E ~ fh_ipolity2 * icrg_qog +ln_cases1 + reg_any_case + weekend + strata(date_firstcase) + as.factor(region) + cluster(ccode) , method = "breslow", data = JPP_covid19_2022_CoxED_E)
summary(covid.cox_E.int)

# E-D-low-QOG:
me_E_int_DlowQOG <- coxed(covid.cox_E.int, method = "npsf",  id = JPP_covid19_2022_CoxED_E$ccode, bootstrap = TRUE, B = 200, maxit = 30,
                          newdata = dplyr::mutate(JPP_covid19_2022_CoxED_E, fh_ipolity2 = JPP_covid19_2022_CoxED_E$lowFPDem, icrg_qog = JPP_covid19_2022_CoxED_E$lowQOG),
                          newdata2 = dplyr::mutate(JPP_covid19_2022_CoxED_E, fh_ipolity2 = JPP_covid19_2022_CoxED_E$highFPDem, icrg_qog = JPP_covid19_2022_CoxED_E$lowQOG))
summary(me_E_int_DlowQOG, stat="mean")

# E-D-high-QOG:
me_E_int_DhighQOG <- coxed(covid.cox_E.int, method = "npsf",  id = JPP_covid19_2022_CoxED_E$ccode, bootstrap = TRUE, B = 200, maxit = 30,
                           newdata = dplyr::mutate(JPP_covid19_2022_CoxED_E, fh_ipolity2 = JPP_covid19_2022_CoxED_E$lowFPDem, icrg_qog = JPP_covid19_2022_CoxED_E$highQOG),
                           newdata2 = dplyr::mutate(JPP_covid19_2022_CoxED_E, fh_ipolity2 = JPP_covid19_2022_CoxED_E$highFPDem, icrg_qog = JPP_covid19_2022_CoxED_E$highQOG))
summary(me_E_int_DhighQOG, stat="mean")

# E-QOG-lowD:
me_E_int_DlowQOG <- coxed(covid.cox_E.int, method = "npsf",  id = JPP_covid19_2022_CoxED_E$ccode, bootstrap = TRUE, B = 200, maxit = 30,
                          newdata = dplyr::mutate(JPP_covid19_2022_CoxED_E, fh_ipolity2 = JPP_covid19_2022_CoxED_E$lowFPDem, icrg_qog = JPP_covid19_2022_CoxED_E$lowQOG),
                          newdata2 = dplyr::mutate(JPP_covid19_2022_CoxED_E, fh_ipolity2 = JPP_covid19_2022_CoxED_E$lowFPDem, icrg_qog = JPP_covid19_2022_CoxED_E$highQOG))
summary(me_E_int_DlowQOG, stat="mean")

# E-QOG-high-D:
me_E_int_DhighQOG <- coxed(covid.cox_E.int, method = "npsf",  id = JPP_covid19_2022_CoxED_E$ccode, bootstrap = TRUE, B = 200, maxit = 30,
                           newdata = dplyr::mutate(JPP_covid19_2022_CoxED_E, fh_ipolity2 = JPP_covid19_2022_CoxED_E$highFPDem, icrg_qog = JPP_covid19_2022_CoxED_E$lowQOG),
                           newdata2 = dplyr::mutate(JPP_covid19_2022_CoxED_E, fh_ipolity2 = JPP_covid19_2022_CoxED_E$highFPDem, icrg_qog = JPP_covid19_2022_CoxED_E$highQOG))
summary(me_E_int_DhighQOG, stat="mean")


## F (E + Additional controls)

JPP_covid19_2022_CoxED_F <- read.csv("JPP_covid19_2022_CoxED_F.csv")

covid.surv_F <- Surv(time = JPP_covid19_2022_CoxED_F$date1, time2 = JPP_covid19_2022_CoxED_F$date2, event = JPP_covid19_2022_CoxED_F$schoolsclosed)
covid.cox_F <- coxph(covid.surv_F ~ fh_ipolity2 + icrg_qog +ln_cases1 + reg_any_case + weekend + ln_gdpcap + wdi_beds_per1000c + wdi_pop14 + wdi_popurb + acd_res_nuclstem + strata(date_firstcase) + as.factor(region) + cluster(ccode) , method = "breslow", data = JPP_covid19_2022_CoxED_F)
summary(covid.cox_F)

ed_F <- coxed(covid.cox_F, method = "npsf",  id = JPP_covid19_2022_CoxED_F$ccode, data=JPP_covid19_2022_CoxED_F)
summary(ed_F, stat="mean")

# F-D:
me_F_D <- coxed(covid.cox_F, method = "npsf",  id = JPP_covid19_2022_CoxED_F$ccode, bootstrap = TRUE, B = 200, maxit = 30,
                newdata = dplyr::mutate(JPP_covid19_2022_CoxED_F, fh_ipolity2 = JPP_covid19_2022_CoxED_F$lowFPDem),
                newdata2 = dplyr::mutate(JPP_covid19_2022_CoxED_F, fh_ipolity2 = JPP_covid19_2022_CoxED_F$highFPDem))
summary(me_F_D, stat="mean")

# F-QOG: 
me_F_QOG <- coxed(covid.cox_F, method = "npsf",  id = JPP_covid19_2022_CoxED_F$ccode, bootstrap = TRUE, B = 200, maxit = 30,
                  newdata = dplyr::mutate(JPP_covid19_2022_CoxED_F, icrg_qog = lowQOG),
                  newdata2 = dplyr::mutate(JPP_covid19_2022_CoxED_F, icrg_qog = highQOG))
summary(me_F_QOG, stat="mean")


## F + Interaction

covid.cox_F.int <- coxph(covid.surv_F ~ fh_ipolity2 * icrg_qog +ln_cases1 + reg_any_case + weekend + ln_gdpcap + wdi_beds_per1000c + wdi_pop14 + wdi_popurb + acd_res_nuclstem + strata(date_firstcase) + as.factor(region) + cluster(ccode) , method = "breslow", data = JPP_covid19_2022_CoxED_F)
summary(covid.cox_F.int)

# F-D-low-QOG:
me_F_int_DlowQOG <- coxed(covid.cox_F.int, method = "npsf",  id = JPP_covid19_2022_CoxED_F$ccode, bootstrap = TRUE, B = 200, maxit = 30,
                          newdata = dplyr::mutate(JPP_covid19_2022_CoxED_F, fh_ipolity2 = JPP_covid19_2022_CoxED_F$lowFPDem, icrg_qog = JPP_covid19_2022_CoxED_F$lowQOG),
                          newdata2 = dplyr::mutate(JPP_covid19_2022_CoxED_F, fh_ipolity2 = JPP_covid19_2022_CoxED_F$highFPDem, icrg_qog = JPP_covid19_2022_CoxED_F$lowQOG))
summary(me_F_int_DlowQOG, stat="mean")

# F-D-high-QOG:
me_F_int_DhighQOG <- coxed(covid.cox_F.int, method = "npsf",  id = JPP_covid19_2022_CoxED_F$ccode, bootstrap = TRUE, B = 200, maxit = 30,
                           newdata = dplyr::mutate(JPP_covid19_2022_CoxED_F, fh_ipolity2 = JPP_covid19_2022_CoxED_F$lowFPDem, icrg_qog = JPP_covid19_2022_CoxED_F$highQOG),
                           newdata2 = dplyr::mutate(JPP_covid19_2022_CoxED_F, fh_ipolity2 = JPP_covid19_2022_CoxED_F$highFPDem, icrg_qog = JPP_covid19_2022_CoxED_F$highQOG))
summary(me_F_int_DhighQOG, stat="mean")

# F-QOG-low-D:
me_F_int_DlowQOG <- coxed(covid.cox_F.int, method = "npsf",  id = JPP_covid19_2022_CoxED_F$ccode, bootstrap = TRUE, B = 200, maxit = 30,
                          newdata = dplyr::mutate(JPP_covid19_2022_CoxED_F, fh_ipolity2 = JPP_covid19_2022_CoxED_F$lowFPDem, icrg_qog = JPP_covid19_2022_CoxED_F$lowQOG),
                          newdata2 = dplyr::mutate(JPP_covid19_2022_CoxED_F, fh_ipolity2 = JPP_covid19_2022_CoxED_F$lowFPDem, icrg_qog = JPP_covid19_2022_CoxED_F$highQOG))
summary(me_F_int_DlowQOG, stat="mean")

# F-QOG-high-D:
me_F_int_DhighQOG <- coxed(covid.cox_F.int, method = "npsf",  id = JPP_covid19_2022_CoxED_F$ccode, bootstrap = TRUE, B = 200, maxit = 30,
                           newdata = dplyr::mutate(JPP_covid19_2022_CoxED_F, fh_ipolity2 = JPP_covid19_2022_CoxED_F$highFPDem, icrg_qog = JPP_covid19_2022_CoxED_F$lowQOG),
                           newdata2 = dplyr::mutate(JPP_covid19_2022_CoxED_F, fh_ipolity2 = JPP_covid19_2022_CoxED_F$highFPDem, icrg_qog = JPP_covid19_2022_CoxED_F$highQOG))
summary(me_F_int_DhighQOG, stat="mean")


##########################################################
# MODELS G & H (Polity: Competition and Participation)
##########################################################

## G

JPP_covid19_2022_CoxED_G <- read.csv("JPP_covid19_2022_CoxED_G.csv")

covid.surv_G <- Surv(time = JPP_covid19_2022_CoxED_G$date1, time2 = JPP_covid19_2022_CoxED_G$date2, event = JPP_covid19_2022_CoxED_G$schoolsclosed)
covid.cox_G <- coxph(covid.surv_G ~ parcomp + parreg + icrg_qog + ln_cases1 + reg_any_case + weekend + strata(date_firstcase) + as.factor(region) + cluster(ccode) , method = "breslow", data = JPP_covid19_2022_CoxED_G)
summary(covid.cox_G)

ed_G <- coxed(covid.cox_G, method = "npsf", id = JPP_covid19_2022_CoxED_G$ccode, data=JPP_covid19_2022_CoxED_G)
summary(ed_G, stat="mean")

# G: Compare Low C and Low P (Congo Brazzaville, 4; 6) to Low C and High P (Belarus,4,8)
me_G_auth_Gl <- coxed(covid.cox_G, method = "npsf",  id = JPP_covid19_2022_CoxED_G$ccode, bootstrap = TRUE, B = 200, maxit = 30, 
                      newdata = dplyr::mutate(JPP_covid19_2022_CoxED_G, parcomp = 4, parreg = 6),
                      newdata2 = dplyr::mutate(JPP_covid19_2022_CoxED_G, parcomp = 4, parreg = 8))
summary(me_G_auth_Gl, stat="mean")

# G: Compare Low C and Low P (Congo Brazzaville, 4; 6)  to High C and High P (Spain, 10, 10)
me_G_auth_Hem <- coxed(covid.cox_G, method = "npsf",  id = JPP_covid19_2022_CoxED_G$ccode, bootstrap = TRUE, B = 200, maxit = 30,
                       newdata = dplyr::mutate(JPP_covid19_2022_CoxED_G, parcomp = 4, parreg = 6),
                       newdata2 = dplyr::mutate(JPP_covid19_2022_CoxED_G, parcomp = JPP_covid19_2022_CoxED_G$highC, parreg = JPP_covid19_2022_CoxED_G$highP))
summary(me_G_auth_Hem, stat="mean")

# G: Compare Low C and High P (Belarus,4,8) to High C and High P (Spain, 10,10) (not reported)
me_G_Gl_Hem <- coxed(covid.cox_G, method = "npsf",  id = JPP_covid19_2022_CoxED_G$ccode, bootstrap = TRUE, B = 200, maxit = 30,
                     newdata = dplyr::mutate(JPP_covid19_2022_CoxED_G, parcomp = 4, parreg = 8),
                     newdata2 = dplyr::mutate(JPP_covid19_2022_CoxED_G, parcomp = 10, parreg = 10))
summary(me_G_Gl_Hem, stat="mean")


## H (G + Additional controls)

JPP_covid19_2022_CoxED_H <- read.csv("JPP_covid19_2022_CoxED_H.csv")

covid.surv_H <- Surv(time = JPP_covid19_2022_CoxED_H$date1, time2 = JPP_covid19_2022_CoxED_H$date2, event = JPP_covid19_2022_CoxED_H$schoolsclosed)
covid.cox_H <- coxph(covid.surv_H ~ parcomp + parreg + icrg_qog + ln_cases1 + reg_any_case + weekend + ln_gdpcap + wdi_beds_per1000c + wdi_pop14 + wdi_popurb + acd_res_nuclstem + strata(date_firstcase) + as.factor(region) + cluster(ccode) , method = "breslow", data = JPP_covid19_2022_CoxED_H)
summary(covid.cox_H)

ed_H <- coxed(covid.cox_H, method = "npsf",  id = JPP_covid19_2022_CoxED_H$ccode, data=JPP_covid19_2022_CoxED_H)
summary(ed_H, stat="mean")

# H: Compare Low C and Low P (Congo Brazzaville, 4; 6) to Low C and High P (Belarus,4,8)
me_H_auth_el <- coxed(covid.cox_H, method = "npsf",  id = JPP_covid19_2022_CoxED_H$ccode, bootstrap = TRUE, B = 200, maxit = 30,
                      newdata = dplyr::mutate(JPP_covid19_2022_CoxED_H, parcomp = 4, parreg = 6),
                      newdata2 = dplyr::mutate(JPP_covid19_2022_CoxED_H, parcomp = 4, parreg = 8))
summary(me_H_auth_el, stat="mean")

# H: Compare Low C and Low P (Congo Brazzaville, 4; 6)  to High C and High P (Spain, 10, 10)
me_H_auth_Hem <- coxed(covid.cox_H, method = "npsf",  id = JPP_covid19_2022_CoxED_H$ccode, bootstrap = TRUE, B = 200, maxit = 30,
                       newdata = dplyr::mutate(JPP_covid19_2022_CoxED_H, parcomp = 4, parreg = 6),
                       newdata2 = dplyr::mutate(JPP_covid19_2022_CoxED_H, parcomp = 10, parreg = 10))
summary(me_H_auth_Hem, stat="mean")

# H: Compare Low C and High P (Belarus,4,8) to High C and High P (Spain, 10,10) (not reported)
me_H_el_Hem <- coxed(covid.cox_H, method = "npsf",  id = JPP_covid19_2022_CoxED_H$ccode, bootstrap = TRUE, B = 200, maxit = 30,
                     newdata = dplyr::mutate(JPP_covid19_2022_CoxED_H, parcomp = 4, parreg = 8),
                     newdata2 = dplyr::mutate(JPP_covid19_2022_CoxED_H, parcomp = 10, parreg = 10))
summary(me_H_el_Hem, stat="mean")
